Lamb waves are guided waves propagating in a solid plate or layer considering stress-free boundary conditions and the propagation of Lamb waves is governed by the Rayleigh–Lamb equation.^[17] By numerically solving the equation, the group velocity dispersion curves for a 2-mm thick aluminum plate with Young’s modulus 69 GPa, density 2700 kg/m³, and Poisson’s ratio 0.33 are plotted in Fig. 1. The multimodal and dispersive propagation behavior of Lamb waves is apparently illustrated. This paper focuses on the within-mode dispersion,^[19] which means different frequency components in a single mode travel at different speeds and result in the spreading of the wave packet as it propagates. Therefore, this characteristic is often referred to as group velocity dispersion.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Group velocity dispersion curves of a 2-mm thick aluminum plate.

In a nonlinear medium, the response of a mono-frequency excitation contains higher harmonic frequencies, which are integer multiples of the excitation frequency. In previous researches, the generation of higher harmonics in micro-damaged waveguides is explained by CAN,^{[10,11,17,20]} where the clapping of interfaces of a micro-structural crack as well as the friction between them lead to asymmetrical dynamic responses and further to higher harmonic generation. CAN is emphasized not only because it is predominant in the case of delamination in composites^[20] but it is also effective in the evaluation of fatigue cracks in metal. A breathing crack with contact analysis is implemented in the following simulations to describe this mechanism. This modeling method is termed “CAN model” in this paper.

Higher harmonic generation is also introduced by the PM model in this paper. Compared with the CAN model, the PM model is the representation of localized material damages at the mesoscopic level. It describes the nonlinear hysteretic behavior of fatigued metallic samples induced by micro-cracks.^[21,22] It is a multiscale model taking the summation of the strain contribution of a large number of elementary hysteretic elements into consideration. The stress–strain relation of each hysteretic element of the PM model considered here is plotted in Fig. 2. Here a pore pressure P, defined as the negative of the local stress, is introduced to describe the hysteretic elements.^[22,23] Each element has two different states, an “open” state and a “closed” state and the transition between these two states is determined by two critical pore pressures P_o and P_c. Assuming the initial state of an element is open, when P is less than P_c and increases, the element represents linear elastic properties with modulus K₁. However, when P leaps over P_c, the strain exhibits a step change with value r₂ and the element changes to a closed state. The element state remains linearly elastic with modulus K₂ afterwards. To the contrary, considering an initially closed element, when P is greater than P_o and decreases, the element represents linear elastic properties with modulus K₂. A step change of strain with value r₁ and a transition to the open state occur when P leaps over P_o. The element state remains open and the element remains linearly elastic with modulus K₁ afterwards. Different hysteretic elements are characterized by different pairs of (P_o, P_c) and the distribution of (P_o, P_c) is described by PM space.^[24] The joint effect of all hysteretic elements over the whole damaged area accounts for the nonlinear dynamics behavior.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. The stress–strain relation of an elementary hysteretic element in the PM model.

Both CAN model and PM model are used in the following simulations. The center of the damaged area is changed to different positions and the location results using the method proposed in this paper are obtained respectively. By comparing the results, two different models of a damaged area responsible for higher harmonic generation are investigated and compared.

The numerical simulations based on nonlinear elastic response and time reversal are performed using ABAQUS/Explicit. As is schematically shown in Fig. 3, two plates with the same size are modeled, dealing with the forward propagation of a mono-frequency excitation signal and the backward propagation of the time-reversed higher harmonics, respectively. The micro-damaged waveguide is an otherwise linear elastic plate with a damaged area (Fig. 3(a)). The linear elastic waveguide is identical with the former one except that it is intact (Fig. 3(b)). The size of these two plates is 1000 × 2 mm and the material parameters are Young’s modulus 69 GPa, density 2700 kg/m³, and Poisson’s ratio 0.33. The CAN model is an ellipse with a major axis of 0.8 mm and a minor axis of 6 nm (Fig. 4(a)). Its major axis is set perpendicular to the surfaces of the plate. The interfaces of the breathing crack are simulated by hard contact with a frictionless model. Moreover, the size of the PM model zone is 2 × 1.5 mm (Fig. 4(b)). The PM material parameters are K₁ = 7.5 × 10¹⁴ Pa, K₂ = 1.0 × 10¹⁵ Pa, r₁ = 0.002, r₂ = 0.001 and consider a uniform distribution of hysteretic elements in PM space.^[25] In order to define PM model material in ABAQUS software, a VUMAT subroutine is implemented to update the stress tensor at every iteration time step.^[22,23,25] The center coordinates of these two models are o_CAN and o_PM, respectively.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The schematic diagram of (a) a micro-damaged plate and (b) an intact plate for the analysis of the forward and backward propagations, respectively.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Two different models of the damaged area: (a) the CAN model and (b) the PM model.

The simulations are performed according to the following steps: (i) S₀-mode Lamb wave is excited in a micro-damaged waveguide and the response signal is recorded at the receiver. (ii) Higher harmonic components induced by the damaged area are extracted and time reversed as the reemitted signal. (iii) Backward propagation of the reemitted signal is performed on an intact waveguide and the signals acquired at the receivers in the inspection region are analyzed.

During the forward propagation, the signal is excited at x = −100 mm and the corresponding response is received at x = 100 mm. The reemission is performed at the same position (x = 100 mm). The reemitted signal is the received signal of the forward propagation processed by PI and time reversal operator. The inspection region of interest is at the center of the upper surface of the intact plate with a size of 100 mm (x = −50 mm to 50 mm). Receivers are set in this region and more dense receivers are placed around the position where the crack is expected for precise location.

Because of the multimodal and dispersive nature of Lamb waves, an accurate determination of higher harmonic components from a time domain signal can be difficult. Compared with short-time Fourier transform whose results strongly depend on setting experience-based parameters, the wavelet transform is a consistent and robust tool once an appropriate mother wavelet is chosen.^[3,26] Therefore, in this paper, the further signal processing for higher harmonic generation is done by the wavelet transform and the Morlet wavelet is employed as the mother wavelet.^[3,17] The peak magnitude of the wavelet coefficient at a given frequency is proportional to the amplitude of the raw signal at this frequency,^[26] hence the higher harmonics can be represented by corresponding wavelet coefficients.

The generation of higher harmonics complicates the nonlinear Lamb wave propagation problem because more frequencies mean more excited Lamb wave modes. However, if an S₀ mode wave is excited at low frequencies in a nonlinear waveguide, it is possible that some of the induced higher harmonics are still S₀ mode and they contain the majority of the harmonic energy. The realization of this S₀–S₀ pair will greatly simplify the problem. Therefore, S₀-mode waves are used in this paper and the excitation of the forward propagation is realized by a hamming windowed tone burst consisting of 5 cycles with the central frequency of 200 kHz.

Furthermore, to ensure the existence of only symmetric modes in the waveguide, symmetric excitation is performed in the simulations of both forward and backward propagations. By simultaneously exciting two sources placed on the double surfaces with the same signal, symmetric mode Lamb waves will be enhanced while antisymmetric waves will be suppressed. In our simulations, the point force method^[27] is used to model the excitation source of Lamb waves. Moreover, the receiver is modeled by a point on the upper surface and the x-direction normal stress σ_xx at this point is recorded as the receiving signal.^[10,23]

With nonlinear analysis as a pretreatment of time reversal, only higher harmonic Lamb waves are expected to be time reversed and reemitted to the waveguide. However, because of the multimodal and dispersive nature of Lamb waves, the extraction of Lamb waves at a certain frequency implicates some difficulties. To overcome this problem, the filtering operation called PI is used. Within a linear medium, the phase inversion of a pulsed excitation signal leads to exact inverted phase response signals. However, it is not the case in a nonlinear medium. It has been proved mathematically that when two phase-inversed excitation signals propagate in a nonlinear medium, in the sum of two corresponding response signals, all of the odd harmonic components including the fundamental component will be cancelled out while all the even harmonic components will be reserved.^[28]

Although PI filtering method has been investigated for bulk waves,^[22] its validation for Lamb wave simulations is essential. In a preliminary simulation, the forward propagation is performed twice using excitation signals with opposite phases (0° and 180°) respectively and the nonlinearity is introduced by CAN model with o_CAN at (0, 0) mm. The response signals are SIG+ and SIG – respectively. Due to the generation of higher harmonics, the sum signal of SIG+ and SIG–, SUM, is a non-zero signal (Fig. 5). Further continuous wavelet transform of these three signals helps identify their characteristics (Fig. 6). As is indicated, all wave packets at odd harmonic frequencies are eliminated while all wave packets at even harmonic frequencies are enhanced in SUM. Among them two wave packets at 400 kHz and 800 kHz, corresponding to the second and fourth harmonic components respectively, are predominate. Furthermore, it is worth noticing that the time-domain distribution of these two wave packets remains unchanged in SIG+, SIG–, and SUM. The simulation results demonstrate that by using PI, it is convenient to extract even harmonic components in nonlinear Lamb waves. Therefore, PI method is employed so that only the relevant information on the localized damage will be extracted and used in our simulations.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Time domain waveforms of (a) SIG+, SIG–, and (b) SUM.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Continuous wavelet transform of (a) SIG+, (b) SIG–, and (c) SUM. All three plots are on the same scale and the color specified by the bar indicates the absolute wavelet coefficient.

PM model is a phenomenological model for micro-damaged materials and has been reported to generate remarkable higher harmonics in bulk wave simulations.^{[22,23,25,29]} The S₀-mode Lamb waves used here resemble the displacement field of the simple axial wave^[30] and show particle displacements mostly in the x direction.^[27] Therefore, higher harmonic generation is expected in nonlinear Lamb wave simulations with the PM model. This is verified by a simulation of forward propagation using PM model with o_PM at (0, 0) mm. Figure 7 shows the wavelet map of the acquired SUM signal under current simulation configuration. As is illustrated, the wave packets at 400 kHz and 800 kHz are clearly visible and still dominant in this case. Comparing Fig. 7 and Fig. 6(c), a similar distribution of these two wave packets in both time and frequency domain can be observed. This means the higher harmonics induced by a localized PM model is comparable to a breathing crack. These observations yield the conclusion that the PM model is an effective model for micro-damages in Lamb wave simulations.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. The wavelet map of SUM signal when a PM model region existing in an otherwise linear elastic plate. The color specified by the bar indicates the absolute wavelet coefficient.

Due to the symmetric excitation performed in both forward and backward propagations, the higher harmonics are constricted to symmetric modes. Furthermore, it is shown by Fig. 1 that at 400 kHz and 800 kHz only the fundamental symmetric mode (S₀) is excited. Therefore, the wave packets at 400 kHz (P₁ packet) and 800 kHz (P₂ packet) are S₀-mode waves during their propagations. Thus the difference between their arrival times at the receiver, as is shown in Fig. 6(c) and Fig. 7, is only attributable to the within-mode dispersion of S₀ mode, which will be compensated by time reversal.^[19]

A preliminary simulation of the backward propagation is performed using time reversed SUM in Fig. 6(c) as the excitation signal. Signals received at x = –50 mm, 0 mm, 50 mm are analyzed and the corresponding wavelet maps are shown in Fig. 8. In addition to P₁ and P₂ which are of interest, higher even harmonics (1200 kHz and 1600 kHz) are clearly recognizable. The group velocities of P₁ and P₂ are estimated respectively, V_P1 = 5178.66 m/s, V_P2 = 4156.28 m/s and marked with red dots in Fig. 1. It is clear that P₁, P₂ are still S₀-mode waves during the backward propagation and the great coincidence of the group velocities indicates that the interaction between these harmonics are negligible. It is also illustrated by Fig. 8 that the only convergence of the harmonic wave packets happens at x = 0 mm, where the damaged area is in the simulation of the forward propagation.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Continuous wavelet transform of signals received at (a) x = –50 mm, (b) x = 0 mm, and (c) x = 50 mm during the backward propagation of time-reversed SUM on an intact plate. All plots are on the same scale and the color specified by the bar indicates the absolute wavelet coefficient.

To locate the damaged area, consider the signal u(x, t) acquired at each receiver in the inspection region during the backward propagation:

The propagation of P₁ and P₂ are investigated to locate the damaged area. Noticing that the S₀–S₀ pair used here does not satisfy the certain conditions of cumulative higher harmonic generation, the higher harmonics are only induced by the localized damage. The damaged area are located by determining the position of E_MAX. To obtain good location results, Δ f = 50 kHz are chosen to calculate e(x, t).

Two groups of simulations are conducted based on the process aforementioned wherein two different models of the damaged area are used, respectively. In the first group of simulations, the damaged area of a plate is modeled by the CAN model. The center coordinate of the crack o_CAN is set to (–20, 0) mm, (–10, 0) mm, (0, 0) mm, (10, 0) mm, (20, 0) mm, respectively. In the second group of simulations, the CAN model is replaced by the PM model. Similarly, the center coordinate of the local defect o_PM is set to (–20, 0) mm, (–10, 0) mm, (0, 0) mm, (10, 0) mm, (20, 0) mm, respectively.

The location results of the CAN model and the PM model are shown in Fig. 9 and Fig. 10 respectively. Relatively accurate locations of the damaged area are obtained when the center of the damaged area is at different positions in both groups of simulations. This is attributed to the significant generation of P₁ and P₂ by using two localized models of the damaged area. Moreover, P₁ and P₂ are extracted and emphasized while their different arrival times at the receiver during the forward propagation are reserved by PI. Additionally, after time reversal, the distribution of harmonic wave packets facilitates the location of the damaged area by just investigating E_MAX, which will locate at the damaged area as it is the source of higher harmonics.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. The location results of the CAN model with its center coordinate oCAN at (a) (–20, 0) mm, (b) (–10, 0) mm, (c) (0, 0) mm, (d) (10, 0) mm, and (e) (20, 0) mm respectively. In each case, the moment (i = 1 – 5) is determined when the maximum of the sum amplitude of P1 and P2, , occurs. The distribution of P1 and P2 at that moment is plotted and the values are normalized by .

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. The location results of the PM model with its center coordinate oPM at (a) (–20, 0) mm, (b) (–10, 0) mm, (c) (0, 0) mm, (d) (10, 0) mm, and (e) (20, 0) mm, respectively. In each case, the moment (j = 1 – 5) is determined when the maximum of the sum amplitude of P1 and P2, , occurs. The distribution of P1 and P2 at that moment is plotted and the values are normalized by .

Table 1 further compares the location results of two different models of the damaged area, CAN model and PM model. It is demonstrated that by combining nonlinear analysis with time reversal and examining the harmonic wave packet distribution, both kinds of defects can be detected and located accurately and the errors of location results are in an acceptable range. Moreover, the numerical results make evident that the CAN model and the PM model are equivalent to describe the nonlinear dynamics behavior of micro-damages in a plate. Therefore, it is reasonable to employ PM model as a phenomenological model of the damaged area in further Lamb wave simulations.

Table 1.

Table 1.

Table 1. The location results of CAN model (1st group) and PM model (2nd group). xd is the center coordinate of the damaged area in the x direction. . xd/mm 1st group 2nd group Result/mm Error/mm Result/mm Error/mm −20 −21.625 −1.625 −20.625 −0.625 −10 −10.875 −0.875 −8.875 1.125 0 −0.857 −0.875 −1.125 −1.125 10 10.125 0.125 10.625 0.625 20 21.375 1.375 19.625 −0.375

Table 1. The location results of CAN model (1st group) and PM model (2nd group). xd is the center coordinate of the damaged area in the x direction. .